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Abstract. 

Co-doped ZnO is the prototypical dilute magnetic oxide showing many of 
the characteristics of ferromagnetism. The microscopic origin of the long range 
order however remains elusive, since the conventional mechanisms for the magnetic 
interaction, such as super-exchange and double exchange, fail either at the fundamental 
or at a quantitative level. Intriguingly, there is a growing evidence that defects both 
in point-like or extended form play a fundamental role in driving the magnetic order. 
Here we explore one of such possibilities by performing ab initio density functional 
theory calculations for the magnetic interaction of Co ions at or near a ZnO {1010} 
surface. We find that extended surface states can hybridize with the e-levels of Co and 
efficiently mediate the magnetic order, although such a mechanism is effective only 
for ions placed in the first few atomic planes near the surface. We also find that the 
magnetic anisotropy changes at the surface from an hard-axis easy-plane to an easy 
axis, with an associated increase of its magnitude. We then conclude that clusters with 
high densities of surfacial Co ions may display blocking temperatures much higher than 
in the bulk. 
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1. Introduction 

Dilute magnetic semiconductors (DMSs) [lj are a new class of materials in which 
transition metal (TM) dopants (usually Co, Fe or Mn) replace native cations in ordinary 
semiconductors. The TMs provide localized spins and in some cases free carriers 
and the material often displays evidence for ferromagnetism at remarkably low TM 
concentrations and relative high temperatures. The expectation around DMSs is that 
they might offer the advantages of semiconductors (i.e. an easy to manipulate electronic 
structure) together with the non-volatility of ferromagnets. In short they may be an 
ideal materials platform for future "spintronics" technologies [2]. 

(Ga,Mn)As is by far the most studied among the DMSs, mostly because of 
its compatibility with present the GaAs/AlAs technology. However, despite more 
than a decade of refinements in the synthetic method and of improvements of the 
annealing treatments to control the various defects concentrations, the (Ga,Mn)As 
Curie temperature still remains far from exceeding room temperature [3], hindering 
the prospect for mainstream applications. These difficulties have stimulated the 
experimental activity towards other host materials and in particular towards oxides 
based DMSs ^ such as ZnO. 

The most commonly used dopant in this case is Co and at present the claims 
for room temperature ferromagnetism in ZnO: Co, following the first exciting report 
from Ueda et al. [5], are numerous (for a review see [6] and references therein). The 
experimental literature however is populated with controversial results and almost 
all possible magnetic phases have been found in samples with the same nominal Co 
concentration and grown under apparently similar conditions. Thus together with 
ferromagnetism, also paramagnetism [TJ and a spin-glass behavior j8] have been reported. 

Already early experiments have pointed to the extreme sensitivity of the magnetic 
properties over the sample preparation, the post-growth treatment and the various 
growth conditions. Deposition protocols at controlled thermodynamical equilibrium, 
such as chemical synthesis [9j [10] and molecular beam epitaxy [TT] generally lead to 
paramagnetic films, while out of equilibrium schemes such as pulsed laser deposition 
by large produce room temperature magnetism [12J. One interesting common feature 
is that in the cases where ferromagnetism is observed the saturation magnetization 
and the remanence are typically small [131 ttS US], while proving a direct correlation 
between magnetism and carriers gave conflicting results [T6| ITT] . More recently, the 
growth of ultra-high quality ZnO: Co films lacking of magnetism even for large n-doping 
conditions [T8l IT9] , indicated clearly that the magnetism must be related to the presence 
of crystallographic defects, either in their point-like or macroscopic form. This line of 
thoughts was used to critically re- analyse a multitude of published data for ZnO:Mn 
and to correlate the appearance of magnetism to a large surface to volume ratio of the 
sample grains [20] . 

With the conventional mechanisms for the magnetic interaction, such as super- and 
double-exchange, failing either at the qualitative or quantitative level [21], most of the 
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theoretical activity in the field has concentrated on first principles methods, in particular 
on density functional theory (DFT). Early calculations for ZnO doped with extremely 
large concentrations of Co (25%) predicted a room temperature ferromagnetic ground 
state in the absence of additional dopants J22J [23] . These concentrations however are 
above the percolation threshold for nearest neighbors magnetic interaction, so that even 
a short range mechanism as super-exchange can, in principle, yield to ferromagnetism. 
Below such a threshold, where essentially all the experiments are conducted, one should 
expect at best only super-paramagnetic clusters j6]. 

Furthermore, it is important to point out that for this rather complicated problem, 
one should take a particular care when choosing the exchange and correlation functional 
to use in DFT. The standard local density and generalized gradient approximations 
(respectively LDA and GGA) severely under-estimate the native ZnO gap and misplace 
in energy the Co d levels. Both effects contribute to predicting spurious long-range 
ferromagnetism. Rigorous corrections based on energy considerations [2U [25j 126] 
improve the description and return a picture of antiferromagnetic interaction among 
nearest neighbors Co ions and little interaction at any other distance [6]. Thus also 
theory points to defects as source of magnetism or at least as a tool to boost either the 
strength or the range of the magnetic interaction. 

For instance uncompensated spins on the surface of secondary CoO phases have 
been proposed as the origin for the experimentally observed magnetic hysteresis 
|27j . This essentially means to attribute the observed ferromagnetism to super- 
paramagnetically blocked clusters. However, this proposal was also rejected at the 
quantitative level by DFT calculations. In fact the calculated order temperatures 
of cubic, wurtzite and zinc-blende CoO polymorphs are all well below room 
temperature (28], with finite size effects not playing any particular role [29J. 




Another possibility is represented by complexes involving Co ions and native ZnO 
intrinsic defects. This was examined in great detailed by Pemmaraju et al [6], who 
concluded that only complexes of O vacancies (Vo) and substitutional Co ions can 
couple magnetically to a medium range if additional n-type doping is present. Other 
defect combinations can provide remarkably large local magnetic interaction but provide 
no long range interaction. A phase diagram was then proposed where the various 
magnetic states could be mapped out onto the relative abundance of Co ions and Co- 
Vo complexes. Although suggestive the presence of Vo in magnetic samples is often 
claimed, it is not clear whether the concentration of Co-Vo complexes needed for room 
temperature ferromagnetism are achievable in reality [30J. 

In this work we look at yet another possibility, i.e. we investigate whether the 
magnetism can arise at a ZnO open surface, namely the non-polar {1010} surface. If 
this is proved true then theory will also converge in establishing a link between the 
macroscopic morphology of samples and their magnetism. 
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2. Calculation Details 

Our calculations are all performed by using DFT. The standard local and semi-local 
exchange and correlation functionals (LDA and GGA) are not appropriate for the 
description of the electronic structure of ZnO and CoO |6J. We have therefore opted 
for an approximated form of the self-interaction correction scheme (ASIC) [31] , which 
has been shown to correctly reproduce the electronic structure of ZnO: Co [6], in good 
agreement with spectroscopical data. The ASIC method is included in a development 
version of the local atomic orbital basis set DFT code Siesta [32] . For all our calculations 
we use the value of a = \ for the ASIC potential scaling parameter. This provides 
a good description of the electronic structure of strongly correlated insulators and 
semiconductors, away from the purely ionic limit |31j. 

The electron density and the Hamiltonian for the valence electrons are described 
by a double-^ polarized (DZP) basis set. Core electrons are replaced by standard norm- 
conserving Troullier-Martins' pseudopotentials. The real space integration grid has a 
spacing equivalent to a 800 Ry plane-wave cut-off, while the reciprocal space sampling 
is performed over a grid with an equivalent real space distance of 15 A. The forces are 
calculated from the LDA Ceperley- Alder [33] functional using the ASIC density [34J, 
and all structures are relaxed to less than 0.01 eV/A. 

3. Un-doped ZnO surface 

The wurtzite structure of bulk ZnO is well described by the ASIC method: both the in- 
plane, a, and the out of plane, c, lattice parameters are slightly contracted with respect 
to experiments. In particular they are calculated to be respectively 3.209 A and 5.180 A 
against the reported values of 3.252 A and 5.313 A. The distortion parameter u after 
relaxation remains essentially identical to its experimental value. 



12 11 10 9 8 7 6 5 4 3 2 1 




Figure 1. Relaxed unit cell used to describe the ZnO {1010} surface: the point of view 
is along the a axis. The numbers label the different atomic layers as used throughout 
in the manuscript. The large purple balls represent Zn, while the red small ones are 
for O. 



The {1010} is the most stable and frequently occurring ZnO [35J surface; it is not 
polar since it contains both Zn and O atoms. The unit cell for the surface is constructed 
from the primitive cell of bulk ZnO and it is illustrated in figure [TJ Due to the non-polar 
nature of the surface we found it unnecessary to passivate the side of the slab opposite 
to the one that gets doped with Co, instead we have chosen to keep the bottom four 
layers of the slab fixed to the bulk crystal structure during relaxation. Furthermore 
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when performing further relaxation we always constrain the cell geometry and allow 
only the atomic positions to relax. A vacuum region of 10 A is sufficient to prevent the 
interaction between slab periodic images, and in fact a further increase of the length 
of the vacuum region produces less than a 1 meV change in the total energy. A slab 
thickness of 17.7 A (12 atomic layers totaling to 24 atoms) is found to be sufficient to 
encompass surface relaxations with a less than 5 meV/A change in the maximum force 
upon introduction of further layers. 

In order to test for surface reconstructions a 2x2 surface cell is created. The 
atoms, excluding those fixed to their ZnO bulk positions, are randomly displaced by a 
displacement of up to 0.5 A. The cell is then allowed to relax, it was found to revert 
back to the unit cell structure (with no reconstruction). Therefore, in the absence of any 
obvious surface reconstruction, we will continue to assume that the primitive surface 
unit cell can be simply constructed from the the primitive ZnO bulk unit cell. 
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Figure 2. Change in the Miilliken charge, Ap, with respect to the ZnO bulk value 
in the first atomic layers of the {1010} surface. The index labeling the position of the 
atomic layers follows the definition of Fig. [U 

In general the formation of the {1010} surface reduces the coordination of both 
O and Zn from 4 to 3 in the top layer, while, of course, the remaining atoms retain 
their bulk coordination number. This results in a modification of both the electronic 
structure and the atomic positions in the vicinity of the surface. In Fig.[2]we present the 
difference, Ap, between the Miilliken charges calculated for the surfaces with respect to 
that of bulk ZnO as a function of the layer position (distance from the surface). We can 
clearly observe that there is a charge re- arrangement in the topmost three atomic layers. 
This is more pronounced at the surface (layer 1), where O increases its Miilliken valence 
charge by 0.03 e (e is the electron charge), mainly due to the enhanced population of 
the 2p orbitals. Of course such an enhancement does not correspond to doping (the 
O p shell is already completely filled) but simply to a re-arrangement of the atomic 
and the overlap Miilliken populations, i.e. to a reduction of the Zn-0 hybridization 
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Figure 3. Electronic structure of the ZnO {1010} surface. In panel (a) we show the 
DOS for the {1010} surface (top panel), bulk ZnO (bottom panel) and bulk ZnO doped 
Co at the Zn site (middle panel). The positive DOSs are for the majority electrons 
while the negative ones are for the minority. The energy region marked in yellow in 
the top panel indicates the energy window used to calculate the local DOS of (b). Note 
the rather localized surface state with amplitude mainly on the O ions. The color code 
for the atoms in panel (b) is the same as that of Fig. [1] 



at the surface. In contrast the Zn Miilliken charge is reduced by 0.06 e, the majority 
of which comes from the 4p orbital. Interestingly the nominal Miilliken population for 
bulk ZnO is already reached after 4 mono-layers for both Zn and O, indicating that the 
perturbation of the electronic structure interests only the surface. 

Further information can be obtained by comparing the density of states (DOS) of 
the surface super cell with that of bulk ZnO [see top panel of Fig. 3(a)| . Together with an 
increase of the Zn d DOS bandwidth the main feature of the surface DOS is the presence 
of a split-off band just at the top of the valence band. This is almost entirely due to O 
p states and it is localized near the surface. In panel Fig. 3(b) we present the isosurface 
corresponding to the local DOS of such a surface state. This corresponds to the charge 
density distribution of energy levels placed within a narrow energy window around the 
valence band top [yellow dashed region in Fig |3(a)j . Clearly the fi gure confirms the 
presence of an O-derived surface state. Finally we note that additional surface states 
appear at the ZnO conduction band bottom and, as a consequence, the ZnO band gap 
is reduced from 3 eV to 2.5 eV. 
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4. Co defects in ZnO 



Cobalt ions enter the ZnO lattice as substitutional defects sitting at the Zn site. Their 
electronic structure and in particular the position of the empty £2 states is still a matter 
of debate [24] , although there is agreement that beyond LDA DFT methods are necessary 
to capture the insulating nature of ZnO:Co. As mentioned previously here we use ASIC 
DFT, which returns an electronic structure in good agreement with UPS data [13J. 
As an example the DOS of a 128 atoms bulk ZnO cell including one single Co ion 



is presented in Fig. 3(a) (middle panel). One can clearly see that the filled minority 



e states lies near the ZnO valence band top, while the first resonant £2 level is well 
within the conduction band [6]. With such an electronic structure n-type doping will 
not open a double exchange channel, while p-type doping is notoriously difficult in ZnO. 
Therefore, as it stands bulk ZnO: Co in absence of other defects cannot sustain a long 
range magnetic interaction. 




Figure 4. Energy of incorporating a substitutional Co defect into the ZnO lattice 
near the {1010}, calculated for a single Co atom in a 2x2 slab. The surface slab for a 
substitutonal Co atom in layer 1 is shown in the insert. 



We now move our attention to investigating the effect of the surface on the 
Co electronic structure and its magnetic interaction. In this case a single Co ion 
is introduced into the slab at different substitutional positions. We use for these 
calculations a 2x2 slab (comprising 96 atoms), which corresponds to a total Zni-^Co^O 
doping of x = 0.21. Note however the the Co concentration in plane is rather high, since 
in a single ZnO atomic layer one out of four Zn is replaced by Co. The cell is always 
relaxed with the exception of the four bottom layers, which are kept at their bulk 
coordinates. In figure H] we plot the energy that a Co ion gains by moving from the bulk 
to the surface. This is calculated as the total energy difference between a given super 
cell with that where Co occupies the middle layer in the slab. A dipole is observed 
for these slab calculations, however this dipole remained relatively constant for all the 
calculations introducing an error of only 5 meV is observed for the worst case, when 
the Co atom was placed in the top layer. Clearly there is a substantial reduction of the 
total energy as Co moves closer to the {1010} surface, indicating that Co segregation 
to the surface is likely if the ions posses enough kinetic energy during the growth. 
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The electronic structure of Co near the {1010} surface is investigated next in 
figure O where we present the DOS of the slab as a function of the Co position. We can 
observe a substantial interaction between the surface states and the Co d levels, both 
at the valence band top and within the conduction band. This is rather strong for the 
two inequivalent positions at the {1010} surface (layer 1 and 2), where a substantial 
broadening of both the e and t 2 peaks is observed. In particular when Co is placed on 
layer 2, which corresponds to the cation site closes to the surface with the full tetrahedral 
coordination, shown in Figure [5j There is a substantial majority d DOS appearing in 
the ZnO gap and in general all the d manifold is pushed rather high in energy, so that 
the native ZnO band-gap is almost entirely filled. This can be understood from the 
distortion of the coordinating oxygen tetrahedron. The Co atom relaxes to a position 
closer to the surface, reducing the bond length to the two surface O to 1.89 A (denoted 
"A" in figure and increasing the bond length to the deepest O to 2.00 A ("B"). 
The bond length of the remaining O ("C") remains unchanged from the bulk value of 
1.96 A. This difference in bond length lifts the degeneracy of the d orbitals and results 
in a smeared DOS. As the Co moves away from the surface the DOS resembles more 
closely a superposition between the DOS of the surface and that of Co doped in the bulk 
of ZnO. Interestingly for all the positions investigated the orbital degeneracy of minority 
e states remain lifted by the interaction with the surface state at the valence band top. 
Again, with the only exception of Co doping in the two first surface layers, such an 
electronic structure is not favorable for sustaining long range magnetic interaction. 




Figure 5. View along b showing the distorted oxygen tetrahedron for a Co atom in 
layer 2. 

5. Surface mediated magnetic coupling 

We next move on to study the magnetic coupling between Co ions placed in the two 
topmost layers close to the {1010} surface. In general we find that the magnetic 
interaction is strongly dependent on the precise mutual position of the defects and 
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Figure 6. DOS and orbital projected DOS (PDOS) for a slab including one single 
Co ion as a function of the distance from the exposed {1010} surface. The number on 
the right-hand side is the layer index, so that 1 corresponds to Co doped in the first 
layer, 2 corresponds to Co doped in the second layer and so on. Note that the surface 
state produce significant changes in the Co DOS only when Co is placed in the first 2 
topmost layers, i.e. at the {1010} surface. 

of the free surface. As such it is necessary to span a large number of possible Co 
positions. In order to perform this analysis we construct from the relaxed super cell 
discussed in section @] a 4x3 surface cell including two substitutional Co ions doped in 
the two topmost layers. In order to make the calculations numerically more feasible our 
constructed cell contains only 6 atomic layers, a number sufficient to keep the forces at 
the surface within 0.01 eV/A. This is equivalent to a overall Co concentration of around 
4% which is a value easily achievable in experiments (again it should be notice that the 
concentration of Co is highly inhomogeneous, i.e. it is very large at the surface). 

The exchange coupling is then calculated with the broken symmetry approach, i.e. 
by comparing the DFT total energy for the same cell, when the magnetic moments 
of the two Co ions are aligned either ferromagnetically, -Efm, or antiferromagnetically, 
Eafi with respect to each other. The energy difference Eaf — -^fm is then a measure of 
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Figure 7. Difference between the total energies of the ferromagnetic (Efm) and 
antiferromagnetic (Eaf) configuration of a cell containing two Co ions as a function 
of the distance between the ions cZco-Co- 

E^ F — E FM is plotted against the distance between the two Co centers, <i Co _ Co - The main 
features are the presence of two distinct regions of strong magnetic coupling, respectively 
ferromagnetic and antiferromagnetic, for short dc -Co, and the almost complete absence 
of long range magnetic interaction. Thus, unless further doped, the {1010} surface does 
not seem to offer a more favorable picture in terms of exchange coupling than the bulk, 
with strong close distance interaction and essentially no long range features. Thus, also 
in this case long range ferromagnetism must be excluded. However, whether or not 
superparamagnetic blocked clusters can form in abundance and mimic a ferromagnetic 
hysteresis signal, remains an open question, whose answer sensitively depends on the 
likely Co concentration achievable on the surface and the Co magnetic anisotropy. This 
will be discussed next and put in prospective in the closing section. 

Going into more details of the magnetic interaction we find that the strong 
ferromagnetic coupling (-Eaf — -^fm > 0) is associated to nearest neighbor Co ions, 
the first being placed right at the surface and the second in the second layer, and 
exchanged coupled through an O also at the surface. In contrast the configurations 
having strong antiferromagnetic interaction {Eaf — -Efm < 0) are characterized by having 
the mediating O atom sitting on the second layer from the surface and the two Co ions 
placed either both on the surface or one on the surface and one in the second layer. This 
suggests that in principle specific pattern of magnetic order should be possible if one is 
able to control the absorption site of the Co ions on the surface. 
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6. Anisotropy 

In general the breaking of the crystal field symmetry at a surface is expected to increase 
the magnetic anisotropy, as for instance is the case of Co on Pt for which the massive 
zero-field splitting of 9.3 meV has been reported [36]. For Co doped in bulk ZnO the 
value for the zero-field splitting is known to be 2.76 cm -1 (0.34 meV) [7], with Co having 
an hard axis along the wurtzite crystallographic c axis and an easy plane (the a-b plane). 
This is consistent with model Hamiltonian calculations for transition metals in trigonally 
distorted tetrahedral coordination [3?]. With such a value in hand and assuming 
a certain degree of frustration we have previously estimated that superparamagnetic 
particles comprising of more than 250 Co ions should be blocked at room temperature 
[6]. The question is now whether magneto-crystal anisotropy at the surface can push 
such size limit towards smaller particles. 

In order to calculate the magneto-crystal anisotropy we use a 2x2 surface cell 
containing 4 atomic layers for a total of 32 atoms, and where one surface Zn ion is 
replaced by Co. An equivalent cell, although with periodic boundary conditions along 
the three dimensions, is used for the bulk calculations. The anisotropy is computed from 
Siesta DFT by using an on-site approximation for the spin-orbit matrix elements [38J. 
This requires a non-collinear calculation, which unfortunately cannot be performed with 
the ASIC method, since at present it is only defined for collinear DFT. Therefore all the 
anisotropy calculations are carried out at the LDA level by using the Ceperley- Alder [33] 
functional. In practice we compute the total energy as a function of the direction of the 
Co magnetic moment with respect to the wurtzite c axis. 




(a) (b) 

Figure 8. (a) Total energy, E, as a function of the angle <f>, between the Co magnetic 
moment and the wurtzite c axis. The calculations are performed for Co in a bulk ZnO 
matrix (ZnO:Co), for wurtzite CoO (WZ CoO) and for a Co ion at the {1010} ZnO 
surface (surface), (b) shows the path taken by the magnetic moment for the surface 
cell. Note </>=0 corresponds to the bulk c direction. 



The results of our calculations are presented in figure 



8(a) 



where we report the total 



Magnetic interaction of Co ions near the {1010} ZnO surface 



12 



energy, E, as a function of the angle, 0, between the magnetization and the wurtzite c 
axis. In the figure we present three cases, namely that of Co in bulk ZnO, of wurtzite 
CoO [29] and the Co doped {1010} surface. For both the bulk ZnO:Co and wurtzite 
CoO, we predict a hard axis along the crystallographic wurtzite c direction, and an easy 
plane perpendicular to this axis, in agreement with experiments. The zero-field splitting 
parameter is estimated by a simply fit to our DFT calculations to be 0.7 meV (assuming 
5 = 3/2 for Co) and it is essentially identical in the two cases. This is approximately 
a factor two larger than what is found in experiment. We are at present uncertain of 
the reasons for such a discrepancy. This may lie in the choice of the LDA functional, 
which usually over-estimates the p-d hybridization, or in the on-site approximation to 
the spin-orbit matrix elements, or in the crude way we use to extract the parameter. 
Even with this uncertainly it is still useful to compare the bulk anisotropy to that of Co 
on the surface. 

We find that the surface changes drastically the nature of the anisotropy. A Co ions 
placed at the under-coordinated position on the surface now has an easy-axis oriented 
at 120° from the wurtzite c cLXIS clS shown in Figure |8(b)| The zero field splitting is 
calculated 1.97 meV, i.e. it is approximately a factor three larger than that in the 
bulk. Assuming that the factor 2 error on our estimate is similar for both bulk and the 
surface, this still gives us a quite remarkably high anisotropy, which pushes the limit for 
superparamagnetic clusters blocked at room temperature to about 80 atoms. This means 
that we need approximately 80 Co atoms, magnetically coupled (i.e. in nearest neighbor 
positions) on the surface, in order to have a cluster displaying magnetic hysteresis at 
room temperature. 



7. Concluding Remarks 

We have studied the electronic structure, magnetic interaction and magnetic anisotropy 
of Co ions at or in the vicinity of the {1010} ZnO surface, in the search for an explanation 
of the often claimed room temperature ferromagnetism in ZnO: Co. Our results give a 
contrasting picture. On the one hand we have found that the magnetic interaction 
is rather strong for atoms that sit on the surface, which is in some circumstances 
ferromagnetic in nature. Moreover there is a substantial energy gain for a single Co 
ion to locate at the surface, so that high superficial Co densities might be expected in 
nanoparticles and granular materials, i.e. in objects with a large surface to volume ratio. 
On the other hand we also find that the magnetic interaction remains extremely short- 
ranged, which means that the surface does not act as an extended defect mediating the 
magnetic coupling. This aspect is not different from what happens in bulk ZnO and 
excludes the presence of long-range room temperature ferromagnetism. Certainly the 
situation might change in the case of polar surfaces or when surface doping is possible. 

There are however a few more facts to consider. Firstly the magnetic anisotropy 
at the {1010} surface is threefold enhanced with respect to its bulk value and it 
changes from hard-axis easy-plane to easy-axis. This means that superparamagnetic 
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clusters formed at the surface will block at much higher temperature than the in bulk. 
Alternatively one can say that for the same blocking temperature, much smaller clusters 
are required at the surface to mimic a ferromagnetic signal. This is indeed an attractive 
prospective which may validate the uncompensated spin model proposed by Dietl et 
al. (27]. Using a rather crude estimate of the size of superparamagnetically blocked 
clusters at room temperature we predict that approximately 80 atoms are required. 
This means that one needs up to 80 atoms connected by a nearest neighbor percolation 
path on the {1010} surface. Considering that the two dimensional percolation threshold 
ranges from 1/2 to about 0.7 depending on the particular surface, we can conclude that 
the majority of the surface Zn sites should be replaced by Co in order to have magnetic 
hysteresis at room temperature. 

This work is sponsored by Science Foundation Ireland and the EU FP7 project 
ATHENA. We thank HPC-Europa and TCHPC for computational support. 
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